Anne Badel & Jacques van Helden
29 mars 2021
Les données, leurs représentations
Comment comparer deux individus
Comment découvrir des “clusters” dans les données ?
Comment déterminer le nombre de groupe optimal ?
Comment comparer deux classifications ?
Les données sont issues de la base Recount2. Nous avons sélectionné l’étude TCGA : The Cancer Genome Atlas, regroupant des données RNA-seq pour plus de 12.000 patients souffrant de différents types de cancer. Nous nous intéressons ici uniquement aux données Breast Invasive Cancer (BIC) concernant le cancer du sein.
Les données ont été préparées pour vous, selon la procédure détaillée au cours sur l’analyse différentielle de données RNA-seq.
Filtrage des gènes à variance nulle et de ceux contenant trop de zéros.
Normalisation (méthode robuste aux outliers)
Analyse différentielle multi-groupes (en utilisant le package Bioconductor edgeR).
Correction des P-valeurs pour tenir compte des tests multiples (nous avons testé ici ~20.000 gènes). Nous estimons le False Discovery Rate (FDR) selon la méthode de Benjamini-Hochberg (fonction R p.adjust(all.pvalues, method="fdr")).
Sélection de gènes différentiellement exprimés sur base d’un seuil \(\alpha = 0.05\) appliqué au FDR.
| X1AB92ADA.637E.4A42.A39A.70CEEEA41AE3 | DA98A67C.F11F.41D3.8223.1161EBFF8B58 | X06CCFD0F.7FB8.471E.B823.C7876582D6FC | A33B2F42.6EC6.4FB2.8BE5.542407A0382E | |
|---|---|---|---|---|
| ENSG00000000003.14 | 16.52775 | 18.84599 | 18.67677 | 17.75458 |
| ENSG00000000419.12 | 17.50349 | 17.01692 | 18.15315 | 19.34170 |
| ENSG00000000457.13 | 17.96040 | 17.29379 | 17.51914 | 18.34040 |
| ENSG00000000460.16 | 16.77151 | 16.30180 | 16.85386 | 17.75886 |
| ENSG00000000938.12 | 15.13011 | 13.57196 | 15.04359 | 15.11192 |
| ENSG00000000971.15 | 17.21898 | 20.05815 | 17.48420 | 18.33249 |
Pour des raisons historiques, en analyse transcriptomique les données sont toujours fournies avec
Cette convention a été établie en 1997, lors des toutes premières publications sur le transcriptome de la levure. Dans ces études, l’objet d’intérêt (l’“individu”) était le gène, et les variables étaient ses mesures d’expression dans les différentes conditions testées.
Pour l’analyse de tissus cancéreux, on considère au contraire que l’“objet” d’intérêt est l’échantillon prélevé sur le patient, et les variables sont les mesures d’expression des différents gènes chez un patient.
Classiquement, en analyse de données, les individus sont les lignes du tableau de données, les colonnes sont les variables.
Ce qui implique de faire attention, et éventuellement de travailler sur la matrice transposée (fonction t() en R) pour utiliser correctement les fonctions classiques.
ENSG00000000003.14 ENSG00000000419.12 ENSG00000000457.13 ENSG00000000460.16 ENSG00000000938.12 ENSG00000000971.15
X1AB92ADA.637E.4A42.A39A.70CEEEA41AE3 16.52775 17.50349 17.96040 16.77151 15.13011 17.21899
DA98A67C.F11F.41D3.8223.1161EBFF8B58 18.84599 17.01692 17.29379 16.30180 13.57196 20.05815
X06CCFD0F.7FB8.471E.B823.C7876582D6FC 18.67677 18.15315 17.51914 16.85386 15.04359 17.48420
A33B2F42.6EC6.4FB2.8BE5.542407A0382E 17.75458 19.34170 18.34039 17.75886 15.11192 18.33249
1 ligne = 1 gène = 1 individu = 1 vecteur
1 colonne = 1 feature = 1 vecteur
l’ensemble des données = 1 data.frame
[1] 819 1000
ENSG00000000003.14 ENSG00000000419.12 ENSG00000000457.13 ENSG00000000460.16 ENSG00000000938.12
X1AB92ADA.637E.4A42.A39A.70CEEEA41AE3 16.52775 17.50349 17.96040 16.77151 15.13011
DA98A67C.F11F.41D3.8223.1161EBFF8B58 18.84599 17.01692 17.29379 16.30180 13.57196
X06CCFD0F.7FB8.471E.B823.C7876582D6FC 18.67677 18.15315 17.51914 16.85386 15.04359
A33B2F42.6EC6.4FB2.8BE5.542407A0382E 17.75458 19.34170 18.34039 17.75886 15.11192
en tenant compte de l’ensemble des individus/ lignes et variables / colonnes = un nuage de points dans un espace à 1000 dimensions
= PAS de représentation possible (pour l’instant)
Nous utiliserons les termes anglais
en français :
On a pas d’information supplémentaire sur nos données, juste le data.frame contenant
Clustering : on cherche à mettre en évidence des groupes (/ des clusters) dans les données
On a une information** supplémentaire sur nos données : on connaît le partitionnement de notre jeu de données
Classification : on cherche un algorithme / un modèle permettant de prédire la classe, le groupe de tout individu dont on connait les caractéristiques
- y a t’il des groupes ? si oui, combien ?
Définition d’une distance : fonction positive de deux variables
Si 1,2,3 seulement: dissimilarité
[1] "MISSING FIGURE"
| Euclidian distance | Correlation coefficient | Correlation distance | |
|---|---|---|---|
| A - B | 4.85 | 0.93 | 0.07 |
| A - C | 5.59 | -0.53 | 1.53 |
| B - C | 1.03 | -0.67 | 1.67 |
dist() avec l’option method = "euclidean", "manhattan", ...| t1 | t2 | t3 | t4 | t5 | |
|---|---|---|---|---|---|
| X | 1.96 | 3.67 | 4.80 | 2.96 | 4.37 |
| Y | 3.83 | 3.29 | 3.28 | 3.59 | 3.34 |
distance euclidienne : dist(mat.xy) = 2.72
distance de manhattan = dist(mat.xy, method = "manhattan") 5.43
1 - cor() avec l’option method = "pearson", "spearman", ...distance de corrélation = 1-cor(t(mat.xy) 1.93
ENSG00000003509.15 ENSG00000001629.9 ENSG00000003096.13 ENSG00000005108.15
ENSG00000001629.9 70
ENSG00000003096.13 70 125
ENSG00000005108.15 58 118 55
ENSG00000006007.11 87 27 141 135
ENSG00000003509.15 ENSG00000001629.9 ENSG00000003096.13 ENSG00000005108.15
ENSG00000001629.9 0.87
ENSG00000003096.13 0.80 0.97
ENSG00000005108.15 0.94 0.74 0.64
ENSG00000006007.11 1.24 0.92 1.28 0.99
X1AB92ADA.637E.4A42.A39A.70CEEEA41AE3 DA98A67C.F11F.41D3.8223.1161EBFF8B58 X06CCFD0F.7FB8.471E.B823.C7876582D6FC A33B2F42.6EC6.4FB2.8BE5.542407A0382E D021A258.8713.4383.9DCA.45E2F54A0411
DA98A67C.F11F.41D3.8223.1161EBFF8B58 0.1296
X06CCFD0F.7FB8.471E.B823.C7876582D6FC 0.0410 0.2803
A33B2F42.6EC6.4FB2.8BE5.542407A0382E 0.0426 0.2437 0.0085
D021A258.8713.4383.9DCA.45E2F54A0411 0.1876 0.4379 0.0899 0.1130
C705FA90.D9AA.4949.BACA.1C022A14CB03 0.0371 0.1702 0.0408 0.0448 0.0774
Revenons à nos données
On peut ensuite essayer de visualiser les données
plot (! ne pas faire si “grosses” données)boxplot (! ne pas faire si “grosses” données)[1] 0
Afin de pouvoir considérer que toutes les variables sont à la même échelle, il est parfois nécessaire de standardiser les données.
soit
soit
classification hiérarchique : mettre en évidence des liens hiérachiques entre les individus
ressemblance entre individus = distance
ressemblance entre groupes d’invidus = critère d’aggrégation
départ : n individus = n clusters distincts
calcul des distances entre tous les individus
regroupement des 2 individus les plus proches => (n-1) clusters
calcul des dissemblances entre chaque groupe obtenu à l’étape \((j-1)\)
regroupement des deux groupes les plus proches => \((n-j)\) clusters
Faire attention au données
Choisir
adaptés à nos données
Les individus dans le plan
=> faire apparaitres des classes / des clusters
construction des centres de gravité des k clusters construits à l’étape \((j-1)\)
\(k\) nouveaux clusters créés à partir des nouveaux centres suivant la même règle qu’à l’étape \(0\)
obtention de la partition \(P_j\)
ENSG00000000003.14 ENSG00000000419.12 ENSG00000000457.13 ENSG00000000460.16 ENSG00000000938.12 ENSG00000000971.15 ENSG00000001036.13 ENSG00000001084.10 ENSG00000001167.14 ENSG00000001460.17 ENSG00000001461.16 ENSG00000001497.16 ENSG00000001561.6 ENSG00000001617.11 ENSG00000001629.9
4 4 4 4 5 1 1 4 1 5 4 1 4 1 1
ENSG00000001630.15 ENSG00000001631.15 ENSG00000002016.17 ENSG00000002330.13 ENSG00000002549.12
1 4 5 4 1
Ces méthodes non supervisées, sont sans a priori sur la structure, le nombre de groupe, des données.
rappel : un cluster est composé
si les individus d’un même cluster sont proches
si les individus de 2 clusters différents sont éloignés => variance inter forte
La coupure de l’arbre à un niveau donné construit une partition. la coupure doit se faire :
après les agrégations correspondant à des valeurs peu élevées de l’indice
avant les agrégations correspondant à des niveaux élevés de l’indice, qui dissocient les groupes bien distincts dans la population.
| k1 | k2 | k3 | |
|---|---|---|---|
| c1 | 57 | 398 | 0 |
| c2 | 0 | 98 | 61 |
| c3 | 293 | 3 | 0 |
| c4 | 0 | 0 | 90 |
| Algorithme | Pros | Cons |
|---|---|---|
| Hiérarchique | L’arbre reflète la nature imbriquée de tous les sous-clusters | Complexité quadratique (mémoire et temps de calcul) \(\rightarrow\) quadruple chaque fois qu’on double le nombre d’individus |
| Permet une visualisation couplée dendrogramme (groupes) + heatmap (profils individuels) | ||
| Choix a posteriori du nombre de clusters | ||
| K-means | Rapide (linéaire en temps), peut traiter des jeux de données énormes (centaines de milliers de pics ChIP-seq) | Positions initiales des centres est aléatoire \(\rightarrow\) résultats changent d’une exécution à l’autre |
| Distance euclidienne (pas appropriée pour transcriptome par exemple) |
Mesure de similarité entre deux clustering
à partir du nombre de fois que les clustering sont d’accord
\[R=\frac{m+s}{t}\]
[1] 0.8153634
\[ \text{ARI} = \frac{\text{RI}-\text{E(RI)}}{\text{Max RI} - \text{E(RI)}}\]
[1] 0.601469
POUR ALLER PLUS LOIN
distance euclidienne ou distance \(L_2\): \(d(x,y)=\sqrt{\sum_i (x_i-y_i)^2}\)
distance de manahattan ou distance \(L_1\): \(d(x,y)=\sum_i |x_i-y_i|\)
distance du maximum ou L-infinis, \(L_\infty\): \(d(x,y)=\max_i |x_i-y_i|\)
distance de Minkowski \(l_p\): \[d(x,y)=\sqrt[p]{\sum_i (|x_i-y_i|^p}\]
distance de Canberra (x et y valeurs positives): \[d(x,y)=\sum_i \frac{x_i-y_i}{x_i+y_i}\]
distance binaire ou distance de Jaccard ou Tanimoto: proportion de propriétés communes
Note : lors du TP, sur les données d’expression RNA-seq, nous utiliserons le coefficient de corrélation de Spearman et la distance dérivée, \(d_c = 1-r\)
Utilisées en bio-informatique:
Distance de Hamming: nombre de remplacements de caractères (substitutions)
Distance de Levenshtein: nombre de substitutions, insertions, deletions entre deux chaînes de caractères
\[d("BONJOUR", "BONSOIR")=2\]
Distance d’alignements: distances de Levenshtein avec poids (par ex. matrices BLOSSUM)
Distances d’arbre (Neighbor Joining)
Distances ultra-métriques (phylogénie UPGMA)
Il existe d’autres mesures de distances, plus ou moins adaptées à chaque problématique :
Jaccard (comparaison d’ensembles): \(J_D = \frac{A \cap B}{A \cup B}\)
Distance du \(\chi^2\) (comparaison de tableau d’effectifs)
Ne sont pas des distances, mais indices de dissimilarité :
Bray-Curtis (en écologie, comparaison d’abondance d’espèces)
Jensen-Shannon (comparaison de distributions) # Distance avec R : indice de Jaccard
ou pour des distances particulières, par exemple l’indice de Jaccard :
| v.a | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
| v.b | 0 | 1 | 0 | 0 | 0 | 1 | 0 |
| v.c | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
v.a v.b
v.b 0.3333333
v.c 0.0000000 0.3333333
par(mfrow = c(1,2))
biplot(prcomp(BIC_expr), las = 1, cex = 0.7,
main = "Données non normalisées")
biplot(prcomp(BIC_expr, scale = TRUE), las = 1, cex = 0.7,
main = "Données normalisées")On considère les données comme des points de \(\mathbb{R}^n\)
\(\mathbb{R}^n\) : espace Euclidien à \(n\) dimensions, où
On considère les données comme des points de \(R^n\) (*)
Sur la base d’une distance
## Plot distances between 3 points in a 2D Euclidian space
plot(x = 0, y = 0, type = "n", xlim = c(0, 5), ylim = c(0, 5),
xlab = "", ylab = "", las = 1,
main = "3 individuals in a 2-D space\nDot plot representation",
panel.first = grid())
points(x = 1, y = 1, col = "blue", pch = 19)
text(x = 1, y = 1, col = "blue", label = "A", pos = 2)
points(x = 2, y = 0, col = "green", pch = 19)
text(x = 2, y = 0, col = "green", label = "B", pos = 4)
points(x = 4, y = 4, col = "red", pch = 19)
text(x = 4, y = 4, col = "red", label = "C", pos = 4)\[D(C_1,C_2) = \min_{i \in C_1, j \in C_2} D(x_i, x_j)\]
\[D(C_1,C_2) = \max_{i \in C_1, j \in C_2} D(x_i, x_j)\]
\[D(C_1,C_2) = \frac{1}{N_1 N_2} \sum_{i \in C_1, j \in C_2} D(x_i, x_j)\]
\(d^2(C_i,C_j) = I_{intra}(C_i \cup C_j)-I_{intra}(C_i)-I_{intra}(C_j)\)
\(D(C_1,C_2) = \sqrt{\frac{N_1N_2}{N_1 + N_2}} \| m_1 -m_2 \|\)
dist() avec l’option method = "euclidean", "manhattan", ...| t1 | t2 | t3 | t4 | t5 | SUM | |
|---|---|---|---|---|---|---|
| X | 4.51 | 2.96 | 4.15 | 3.43 | 1.12 | 16.19 |
| Y | 3.33 | 2.81 | 2.09 | 3.26 | 2.65 | 14.14 |
| abs(Y - X) | 1.18 | 0.16 | 2.06 | 0.18 | 1.53 | 5.11 |
| (Y - X)^2 | 1.39 | 0.02 | 4.26 | 0.03 | 2.34 | 8.05 |
| Eucl | 1.18 | 0.16 | 2.06 | 0.18 | 1.53 | 2.84 |
distance euclidienne : 2.84
distance de manhattan = 5.11
1 - cor() avec l’option method = "pearson", "spearman", ...distance de corrélation = 0.84
… c’est à dire aux options des fonctions dist() et hclust()
par(mfrow = c(2, 1))
plot(BIC.hclust, hang = -1, cex = 0.5, main = "Données brutes")
plot(BIC.scale.hclust, hang = -1, cex = 0.5, main = "Normalisées")| # R environment used for this analysis |
r ## Print the complete list of libraries + versions used in this session sessionInfo() |
| ``` R version 4.0.2 (2020-06-22) Platform: x86_64-apple-darwin17.0 (64-bit) Running under: macOS Mojave 10.14.6 |
| Matrix products: default BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib |
| locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 |
| attached base packages: [1] stats graphics grDevices utils datasets methods base |
| other attached packages: [1] tinytex_0.27 pheatmap_1.0.12 vegan_2.5-6 lattice_0.20-41 permute_0.9-5 RColorBrewer_1.1-2 aricode_1.0.0 FactoMineR_2.3 knitr_1.30 |
| loaded via a namespace (and not attached): [1] Rcpp_1.0.5 highr_0.8 pillar_1.4.6 compiler_4.0.2 tools_4.0.2 digest_0.6.27 nlme_3.1-150 evaluate_0.14 lifecycle_0.2.0 tibble_3.0.4 gtable_0.3.0 mgcv_1.8-33 pkgconfig_2.0.3 rlang_0.4.8 [15] Matrix_1.2-18 parallel_4.0.2 ggrepel_0.8.2 yaml_2.2.1 xfun_0.19 stringr_1.4.0 dplyr_1.0.2 cluster_2.1.0 generics_0.1.0 vctrs_0.3.4 flashClust_1.01-2 grid_4.0.2 tidyselect_1.1.0 scatterplot3d_0.3-41 [29] glue_1.4.2 R6_2.5.0 rmarkdown_2.5 farver_2.0.3 ggplot2_3.3.2 purrr_0.3.4 magrittr_1.5 splines_4.0.2 scales_1.1.1 ellipsis_0.3.1 htmltools_0.5.0 leaps_3.1 MASS_7.3-53 colorspace_2.0-0 [43] labeling_0.4.2 stringi_1.5.3 munsell_0.5.0 crayon_1.3.4 ``` |
Contact: anne.badel@univ-paris-diderot.fr
Comment comparer des vecteurs-individus ?